﻿// #define DEBUG_OUTPUT

using PengSW.RuntimeLog;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace PengSW.Math
{
    public static class LinearEquations
    {
        public static bool IsZero(double x)
        {
            return x >= -double.Epsilon && x <= double.Epsilon;
        }

        /// <summary>
        /// 矩阵行调整
        ///     将源行SourceRow乘以因子Delta后，加到目标行DestRow上。
        /// </summary>
        /// <param name="M">矩阵</param>
        /// <param name="SourceRow">源行标</param>
        /// <param name="DestRow">目标行标</param>
        /// <param name="Delta">因子</param>
        private static void MatrixRowAdjust(double [,] M, int SourceRow, int DestRow, double Delta)
        {
            Parallel.For(0, M.GetLength(1), i => M[DestRow, i] += M[SourceRow, i] * Delta);
        }

        /// <summary>
        /// 在终端输出指定矩阵，用于调试输出。
        /// </summary>
        /// <param name="M">目标矩阵</param>
        public static void PrintMatrix(double[,] M)
        {
            int n = M.GetLength(0);
            int m = M.GetLength(1);
            StringBuilder aMatrixString = new StringBuilder();
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j < m; j++)
                {
                    aMatrixString.Append(string.Format("{0}, ", M[i, j]));
                }
                aMatrixString.AppendLine();
            }
            RL.WriteLog(aMatrixString.ToString());
        }

        public static double[,] Elimination(double[,] aExpandMatrix)
        {
            int rowcount = aExpandMatrix.GetLength(0);
            int colcount = aExpandMatrix.GetLength(1);
            if (rowcount >= colcount) throw new System.ApplicationException("增广矩阵的列数必须大于行数！");
#if DEBUG_OUTPUT
            PrintMatrix(aExpandMatrix);
#endif
            // 逐一对对角线元素[i, i]进行消元
            for (int i = 0; i < rowcount; i++)
            {
                // 如果对角线元素为零，则在以下行的这一列中找出最大元素整行加到i行上
                if (IsZero(aExpandMatrix[i, i]))
                {
                    // 找此列中的最大值，加到对角线元素上
                    double aMaxAbsValue = 0;
                    int k = i;
                    for (int j = i + 1; j < rowcount; j++)
                    {
                        double aAbsValue = System.Math.Abs(aExpandMatrix[j, i]);
                        if (aAbsValue > aMaxAbsValue)
                        {
                            aMaxAbsValue = aAbsValue;
                            k = j;
                        }
                    }
                    if (IsZero(aMaxAbsValue)) throw new System.ApplicationException("增广矩阵的秩为0，方程组无解！");
                    Parallel.For(0, colcount, col => aExpandMatrix[i, col] += aExpandMatrix[k, col]);
#if DEBUG_OUTPUT
                    RL.WriteLog(string.Format("第{0}列最大值为{1}行的{2}.", i, k, aMaxAbsValue));
                    RL.WriteLog(string.Format("将第{0}行加到{1}行上：", k, i));
                    PrintMatrix(aExpandMatrix);
#endif
                }
                // 用第i行对其他各行进行消元
                double aDelta;
                for (int j = 0; j < rowcount; j++ )
                {
                    if (i == j) continue;
                    aDelta = -aExpandMatrix[j, i] / aExpandMatrix[i, i];
                    Parallel.For(0, colcount, col => aExpandMatrix[j, col] += aExpandMatrix[i, col] * aDelta);
#if DEBUG_OUTPUT
                    RL.WriteLog(string.Format("用第{0}行对第{1}行消元：", i, j));
                    PrintMatrix(aExpandMatrix);
#endif
                }
                aDelta = aExpandMatrix[i, i];
                Parallel.For(i, colcount, col => aExpandMatrix[i, col] /= aDelta);
#if DEBUG_OUTPUT
                RL.WriteLog(string.Format("将[{0},{0}]置1：", i, i));
                PrintMatrix(aExpandMatrix);
#endif
            }
            // 取出结果矩阵
            double[,] aResultMatrix = new double[rowcount, colcount - rowcount];
            for (int i = 0; i < rowcount; i++) for (int j = 0; j < colcount - rowcount; j++) aResultMatrix[i, j] = aExpandMatrix[i, j + rowcount];
            return aResultMatrix;
        }
    }
}
